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We investigate the zero-temperature phase diagram of interacting Bose gases in the presence of a 
simple cubic optical lattice, going beyond the regime where the mapping to the single-band Bose- 
Hubbard model is reliable. Our computational approach is a new hybrid quantum Monte Carlo 
method which combines algorithms used to simulate homogeneous quantum fluids in continuous 
space with those used for discrete lattice models of strongly correlated systems. We determine the 
critical interaction strength and optical lattice intensity where the superfluid-to- insulator transition 
takes place, considering also the regime of shallow optical lattices and strong interatomic interactions. 
The implications of our findings for the supersolid state of matter are discussed. 
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In quantum many-body systems, the interplay between 
intcrparticlc interactions and commensurability effects 
in external periodic potentials can induce a transition 
from a conducting (or superconducting) to an insulating 
phase, where particles are localized around the minima of 
the periodic potential. This phenomenon (known as the 
Mott transition) is responsible for the insulating behavior 
of several transition metal compounds [l[ . Recently, the 
Mott transition has been induced in a controlled way in 
both bosonic 0] and fermionic [HIH atomic gases trapped 
in optical lattices, which are being used as quantum em- 
ulators to explore the intriguing properties of strongly 
correlated materials Q. 

Previous theoretical studies of cold atomic gases in op- 
tical lattices are based on single-band lattice models, such 
as the Bose- or Fcrmi-Hubbard models Q . These theories 
arc limited to the regime of large optical lattice intensity - 
since excited Bloch bands are neglected - and weak in- 
teractions, since the true interatomic potential is replaced 
by a nonregularized 5 function, which corresponds to 
the first Born approximation for the scattering problem. 
The role of higher Bloch bands has been considered in 
Refs. lZ|,|8| (by using the nonregularized pscudopotential), 
resulting in effective single-band models with density- 
dependent Hubbard parameters or many-body interac- 
tions. The description of interactions beyond the first 
Born approximation is a challenge, since more accurate 
pseudopotentials such as the regularized 5 function in- 
duce virtual transitions to an extremely large number of 
Bloch bands 0. The validity of the single-band approx- 
imation has been ques tioned on an even more profound 
level by Anderson |ToL [Tlj. in particular, concerning the 
existence of a bosonic Mott insulator. One central issue is 
the presence of nodes in the ground-state wave function 
of the bosonic Mott phase in the single-band Hubbard 
model, since this is built on the basis of orthogonal -and, 
hence, nonpositive definite - Wannier functions. Ander- 
son further argues that when all Bloch bands are included 
the superfluid component is always finite, implying that 




5-1(T 0.01 



0.02 0.04 0.08 0.16 0.32 0.64 



a s /d 



FIG. I: (color online). Zero-temperature phase diagram as 
a function of the optical lattice intensity Vo/E T and the s- 
wave scattering length a s /d at unit filling. E r is the recoil 
energy and d is the optical lattice periodicity. The red points 
indicate the transition from the superfluid (green region) to 
the insulating phase (blue region) for hard-sphere particles. 
The solid gray curve is the phase boundary of the single-band 
Bose-Hubbard model [2^| (for a nonregularized 5-function in- 
teraction [(|). The dashed yellow curve is the multiorbital 
model of Ref. (obtained within the mean-field approxi- 
mation). The black diamonds are the experimental results 
of Ref. JT3]. The solid red curve is a guide to the eye, and 
the dashed red curve is an extrapolation towards the freezing 
point without an optical lattice (red vertical arrow). See the 
text for other scenarios that are possible in this regime. 



the solid phase is a supersolid. 

In this Letter, we map out the zero-temperature and 
continuous-space phase diagram of a Bose gas with short- 
range interactions in a simple cubic optical lattice at a 
density of one particle per unit cell (see Fig.Q]), by using 
a new numerically exact hybrid quantum Monte Carlo 
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technique which permits us to simulate the ground state 
of interacting Bose gases in arbitrarily weak or strong pe- 
riodic potentials. Our main goal is to uncover the physics 
of the Mott transition in this continuous-space system, 
going beyond approximate discrete lattice models. We 
find a superfluid-to-insulator quantum phase transition, 
calculate the critical optical lattice intensity as a function 
of the s-wave scattering length, and compare our results 
to an experiment performed with an ultracold gas of Cs 
atoms jl2l ]. finding good agreement (see Fig. [1]). 

Bosons in an optical lattice are described by the 
continuous-space Hamiltonian 
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[V? + V(rO]+$>(ry), (1) 
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where m is the atomic mass, the indices i and j label the 
N particles at positions r^, and ry = — rj. The op- 
tical lattice potential is V (r) = Vq ^2 a=x y z sin 2 (wa/d), exp 
where Vq is the intensity of the laser beams that create 
the lattice and d is half the wave-length. 

The third term in Eq. (fT]) describes the two-body in- 
teractions. Different models for the interparticle inter- 
actions can be considered, including potentials with a 
hard core or with long-range interactions. One spe- 
cific potential we use is the hard-sphere (HS) potential 



(r) = oo, if r < a s and otherwise. The diameter 
of the sphere is equal to the s-wave scattering length a s , 
which is the parameter that characterizes the strength 
of interactions in dilute gases at low temperature. To 
quantify the role of other details of the interparticle po- 
tential, which become relevant at high scattering energies 
and at short interparticle distances, we also employ dif- 
ferent models with the same s-wave scattering length, 
the negative-power potentials (NP) i> NP (r) = c/r p , with 
c > and the integer p > 3 [l3( ■ We do not consider at- 
tractive potentials supporting bound states, such as those 
used to model Feshbach resonances [3] , since in this case 
the possible metastable gaslike phase is not the ground 
state. 

To map out the phase diagram, we employ a new algo- 
rithm combining the stochastic integration of variational 
trial wavefunctions defined on a discrete lattice with an 
imaginary-time propagation in continuous space. The 
ground-state properties of a homogeneous quantum fluid 
[V(r) = 0] can be obtained by using the ground-state 
path- integral Monte Carlo method [15] . where the ex- 
pectation value of a generic operator O is obtained as 

6 = J dRdR/^ (R) e-^Oe-^ T (RQ 
/ dRdR'** (R) e-P^T (R') 

where R = (n, . . . ,rjv). For a sufficiently long imagi- 
nary time j3, this expression gives the exact ground-state 
result provided that the trial wave function "fx (R) has a 
finite overlap with the ground state. The imaginary-time 

evolution operator exp I —/3/2H) can be represented in 



the form of a discrctized path integral using the Trotter 
break up, and the resulting multivariates integral can be 
evaluated by using standard Monte Carlo sampling tech- 
niques [l6l. |22|] . However, the simulation of the strongly 
correlated regime in a periodic potential with the ground- 
state path-integral Monte Carlo method is unfeasible 
with standard trial wave functions. We thus combine the 
ground-state path-integral Monte Carlo algorithm with 
the stochastic sampling of variational wave functions de- 
fined on a discrete lattice. The lattice wave function is 
characterized by the occupation numbers of a set of or- 
bitals localized at the nodes of the lattice. Determining 
the expectation value ([2]) with such trial wave functions 
requires the sampling of those quantum numbers together 
with the continuous-space coordinates R, as explained 
below. 

As a lattice wave function we use a bosonic Gutzwiller 
ansatz with contact interaction [It], [3 <j> (li, • • • , Ijv) = 
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where 5^ is the Kronecker delta, 

7 is a variational parameter, and 1 refers to one of 
the sites of a simple cubic lattice. For the local- 
ized orbitals, we consider a (positive definite) Gaus- 
sian approximation of the Wannicr functions wi(r) = 
A/" 1 / 2 ^ k exp (— zk ■ 1) &k(r), where M is the number of 
lattice sites and 6k (r) is the Bloch state of the lowest 
band with quasimomentum k. For a given lattice con- 
figuration L = . . . , ljy), the coordinate-space wave 
function corresponding to <f> (L) is ipj, (ri, . . . , rjy) = 
n»=i' lo l* ( r «')i so that the combined trial wave function 
is *t (R) = J2l, ipL (R), where the sum goes over 
all lattice configurations L. To improve beyond the 
contact-interaction ansatz we introduce a Jastrow term 
X (R) = Ili<j / ( r jj )i giving a final trial wave function 
f T (R) = x(RVEi >(L) V> L (R)- The Jastrow corre- 
lation function [ljj / (r) is set equal to the free-space 
solution of the two-body problem up to a cutoff distance 
r c and equal to 1 for |r| > r c , as in Ref. [2l|. With this 
choice, the short-range correlations due to the hard core 
of the interatomic interaction v (r) are treated exactly in 
the trial state ^t, while the long-range correlations can 
be efficiently reproduced by the imaginary-time propaga- 
tion. The time evolution also includes the effect of the 
higher Bloch bands. 

To evaluate the right-hand side of Eq. ([2]), we need 
to perform a stochastic sampling over the lattice config- 
urations ^ L and the particle coordinates J dR. This 
requires a random walk in the combined configuration 
space {L, R}, which we perform by using local Monte 
Carlo updates on the lattice index 1^ and the space coor- 
dinate Ti of a given particle i. 

In the present problem, the bosonic statistics plays a 
central role. The trial state is by construction sym- 
metric under particle exchange, and the imaginary-time 
propagation is symmetrized by performing particle per- 
mutations using the worm algorithm [22| . These permu- 
tations can efficiently introduce more exotic short-range 
correlations which are not explicitly included in the trial 
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FIG. 2: (color online). Main panel: Condensate fraction vs 
optical lattice intensity for HS with a a — 0.3d (blue bullets). 
The black solid line is a linear fit used to extract the critical 
point Vo ■ The yellow square at Vo = is the diffusion Monte 
Carlo result of Ref. [2l| . Inset: Finite-size scaling analysis of 
no; L is the simulation-box size. The black solid line is the 
scaling behavior expected at the critical point: no oc L~ 2 . 



wave function * T [HHl- 

For some specific values of 7, the trial state \&t re- 
produces wave functions that have previously been used. 
When 7 — > 0, converges to the "periodic+Jastrow" 
wave-function * T (R) ^ X (R) IL b *=o 0] (ignor- 
ing an irrelevant normalization coefficient). This wave 
function accurately describes the ground state of a su- 
pcrfluid Bose gas in a periodic potential. In the oppo- 
site limit 7 — >• co. one obtains a symmetrized state with 
one particle per site [1H (we consider the case TV = M) 

*t(R) ^ x(R)Ep(L)ll l w p(i I ) ( r i)> where the sum 
runs over the iV! permutations P(L) = (p(li), . . . , p(l/v)) 
of N particles in M = N lattice sites. In numerical stud- 
ies of bulk helium-4, a state analogous to this was intro- 
duced in order to describe a potential supersolid ground 
state [ID, [23], where superfluidity coexists with diago- 
nal long-range order (in this context it is referred to as 
"permanent + Jastrow" , and the Wannier functions are 
replaced by effective localized orbitals). However, it has 
never been implemented in other exact quantum Monte 
Carlo methods such as the diffusion Monte Carlo sim- 
ulations because of the difficulty of performing random 
walks in permutation space. In the context of the single- 
band Bose-Hubbard model, the Gutzwiller ansatz (L) 
predicts, in the limit 7 — > 00, a quantum phase transition 
from a superfluid to an insulating solid [13, 13 • This 
transition is also seen by exact lattice quantum Monte 
Carlo simulations [28j . Here we study this transition in 
the full Hilbcrt space of the continuum model and inves- 
tigate deviations from the single-band lattice model. 

We identify the superfluid phase by measur- 
ing the condensate fraction [20], which is effi- 
ciently evaluated in the worm algorithm as no/n = 
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FIG. 3: (color online). Main panel: Condensate fraction 
vs s-wave scattering length <z s for different model potentials. 
The optical lattice intensity is Vo = 3-Er. Inset: Condensate 
fraction as a function of the variational parameter 7. Here 
Vo — 3.5-Br and a s — 0.3d. The curves corresponding to dif- 
ferent imaginary times p cross at the optimal 7. 



TV" 1 Jdrlim| r _ r /|_ >00 (>l' t (r') * (r)}, where #t (#) is the 
creation (annihilation) field operator. The results for the 
HS potential with a s = 0.3d are shown in Fig. [21 In 
the homogeneous limit Vo = 0, we find agreement with 
the diffusion Monte Carlo prediction of Ref. [23 (yellow 
square). By increasing Vo, the condensate fraction de- 
creases monotonically. Close to a critical optical lat- 
tice intensity (call it V C ), uq approaches zero linearly, 
consistently with the scaling law of the order param- 
eter in the Bosc-Hubbard model (^(r)) oc \fno~Jn oc 

[[Vo — Vn C )/Vn C ] , where /3 = 1/2 (plus logarithmic cor- 
rection) (29[. With a linear extrapolation to no = 
(black solid line in Fig. [2]), we obtain the critical value 
V§/E T = 3.95(10), where E T = h 2 ir 2 /2md 2 is the recoil 
energy. Different extrapolation procedures yield com- 
patible results. In order to better identify the quantum 
phase transition, we perform a finite-size scaling analy- 
sis. In the inset in Fig. [21 we show no/n as a function of 
the simulation-box size L (we employ periodic boundary 
conditions) for different values of Vo . While slightly away 
from V C the condensate fraction rapidly saturates to the 
thermodynamic limit, close to the critical point the scal- 
ing of no/n approaches the power law no oc L~ 2 , which 
is indeed the scaling expected for the order parameter 
squared [29j . 

In the vicinity of the critical point V C , the imaginary 
time /3 required to converge to the ground state scales 
with the size of the system. Hence, an accurate opti- 
mization of "fx is important to avoid a possible bias due 
to the variational ansatz. Rather than by minimizing 
the ground-state energy, this optimization is achieved 
more conveniently with a finite imaginary-time scaling 
analysis, as shown in the inset in Fig. [31 The values of 
no/n obtained with a small imaginary time /3 show a 
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clear dependence on the variational parameter 7. This 
dependence vanishes for larger /3. The crossing point of 
the curves of no /n vs 7 corresponding to different values 
of /? provides the optimal value of 7 which removes the 
residual hnite imaginary-time effects. 

In the phase diagram of Fig. [TJ we show V Q C as a func- 
tion of a s . In the region Vq < V ° , the condensate fraction 
is finite and the system is supcrfluid. The region where 
?io = is instead identified as an insulating phase. The 
phase boundary obtained with the HS potential (red cir- 
cles) is compared to that of the Bose- Hubbard model (28j 
(lower solid gray curve), based on the mapping to single- 
band models described in Ref. Jf|. For a s < 0.3d and 
Vq > 4E r , the two models follow the same trend, with 
small discrepancies due to the effect of higher Bloch 
bands and to the different models for the interatomic po- 
tential (HS vs nonregularized delta). The experimental 
data [12( (black diamonds) are consistent with both the 
continuum and lattice model but are in better agreement 
with our continuum model. The yellow dashed curve 
in Fig. [T] corresponds to a Bose-Hubbard model with 
occupation-numbcr-dcpendcnt parameters due to multi- 
orbital effects, solved within the mean-field approxima- 
tion (from Ref. Q ) , and clearly simulations going beyond 
mean-field will be important to compare this model to 
our continuum results. 

In the regime a s > 0.3d, the HS result bends down, 
as one approaches the freezing density of hard spheres in 
free space V(r) = [3(| (red vertical arrow). This re- 
gion of the phase diagram could host more exotic phases, 
like a quantum-glassy phase eventually coexisting with 
superfluidity, due to the competition between the spatial 
symmetry of the optical lattice (simple cubic) and that 
of the hard-sphere crystal in free space (face-centered cu- 
bic). We leave the investigation of these issues to future 
work. 



In this part of the phase diagram, universality in terms of 
a s is lost and also other details of the interatomic poten- 
tial become important. We quantify these nonuniversal 
effects in Fig. [3] (main panel) by comparing the results for 
no/n vs a s /d (at fixed Vq = 3E r ) obtained with different 
model potentials. The data corresponding to HS (black 
squares) and to NP with n = 9 (blue circles) agree within 
10% up to a s = 0.3d, indicating a large degree of univer- 
sality As expected, the data for NP with n — 6 (cyan 
triangle) show larger deviations. Notice that for n < 6 
the effective range diverges [3l|. In the limit a s — > 0, 
no/n converges to no/n = 0.833, which is the square of 
the overlap between the condensate wave function (the 
Bloch state with k = 0) and the plane wave with zero 
momentum. 

Our results for the Mott transition in a continuous- 
space system of bosons with short-ranged interactions 
can be useful to extend experiments performed with ul- 
tracold gases to the regime of strong interactions and 
weak optical lattices, that is beyond the physics of the 
Hubbard model. While the existence of a quantum 
phase transition in the single-band model was rigorously 
proven [32| . our findings indicate that quantum critical 
behavior persists when all Bloch bands are included, in 
contrast to the arguments of Anderson [h], [ll| ■ 

The hybrid lattice-continuum quantum Monte Carlo 
method presented in this work can also be extended 
to fermionic models and be combined with fixed-node 
and transient-estimate techniques. This will provide a 
new approach for the simulation of strongly correlated 
fermionic systems in weak optical lattices. 
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